		  
#include "stdafx.h"
using namespace std;

     Mod_wc::Mod_wc(Mod_in *iinp):Mod_out_handling(iinp){;}
    Mod_wc::~Mod_wc(){;}
   int Mod_wc::compute(time_intervall &ti){}
    
Mod_wc_fc::Mod_wc_fc(map <string,int> &iIO_options,
                  Mod_in *iinp,
                  Mod_PDF *imod_pdf,
                  Mod_ETP *imod_etp,
                  Mod_Soiltemp *imod_soiltemp,
                  Mod_Surface_water *imod_surface_water,
                  Mod_Crop *imod_crop):Mod_wc(iinp)
  {
		//1. switches setzen
    IO_options=iIO_options;
    //2. Input_objekt setzen
  
    //3. Parameter definieren
    epauft=0.01;// Aufteilung der ETP auf einzelene layer
    
    
    //4. Modelle verkn�pfen
    mod_pdf=imod_pdf;
    mod_etp=imod_etp;
    mod_crop=imod_crop;
    mod_soiltemp=imod_soiltemp;
    mod_surface_water=imod_surface_water;

    
    //5. Output_kategorien f�r diese Variante setzen!!

    this->define_field_of_states("wc",inp->horizons->size());
    this->define_field_of_fluxes("q",inp->horizons->size()-1);
    
    timepointmap <double> nmap (inp->tempgrid);
    this->out->states["Potential_Infiltration"]=nmap;
    //this->out->states["max_EP"]
   
    
  
  }
int Mod_wc_fc::compute(time_intervall &ti)
{
	compute_one_waha_cycle(ti);
	return 0;
}
int Mod_wc_fc::compute_one_waha_cycle(time_intervall &ti)
{
    int zn;
    int x1;
    double real_infiltration;
    timestamp start=ti.get_start();
    timestamp end=ti.get_end();
		for (zn=0; zn<this->out->field_of_states["wc"].size(); zn++) // anfang reinkopierter Teil
		{
			
			x1=(zn==0)+(zn==this->out->field_of_states["wc"].size()-1)*2;
			switch (x1)
			{
				case 1:// oberstes Kompartiment
					
					real_infiltration=wcint (zn,ti);			   // [cm/timeintercall]
					//flux from surface to first layer
					this->set_push_back_fluxes("q",zn,ti,real_infiltration*10/ti.get_dur());		//[mm]
					//q[zt][zn]=max_float(net_net,(wcs[zt][zn]-wc[zt-1][zn])*dx);//Bodenwasserstrom nach unten ist positv
					this->set_states("wc",zn,end,
										this->get_states("wcs",zn,end)
										-this->guelpprog(zn,ti)
										*(this->mod_soiltemp->get_states("t",end)>0)
										/this->get_dx(zn));
										
										
					//wc [zt] [zn]=  wcs [zt] [zn]-(guelpprog(x,zt, zn+1))*(sote[zt][zn]>0)/dx;
					
					//printf("%i  %i  %f\n", zt, zn);
				break;
				case 2:// unterstes Kompartiment

					wcint (zn,ti);	// calculates (wcs,zn,end)
					this->set_push_back_fluxes("q",zn,ti,							 
												Basics_math::max_(0.0,
												(this->get_states("wcs",zn,end)
												-this->get_states("wc",zn,start))
												*this->get_dx(zn)*10/ti.get_dur()));			//flux from, zn-1 to zn   [ mm/d]
					//q [zt] [zn-1]= max_float(0,(wcs[zt][zn]-wc[zt-1][zn])*dx);
					this->set_states("wc",zn,end,
												this->get_states("wcs",zn,end)
												-guelpprog(zn+1,ti)/this->get_dx(zn)
												*(this->mod_soiltemp->get_states("t",zn,end)>0));
					//wc [zt] [zn]= wcs [zt] [zn]-guelpprog(x,zt, zn+1)*(sote[zt][zn]>0)/dx;
					this->set_push_back_fluxes("q",zn+1,ti,
												Basics_math::max_(0.0,
												this->get_states("wcs",zn,end)
												-this->get_states("wc",zn,end))
												*this->get_dx(zn)*10/ti.get_dur());		 //flux from, zn to zn+1   [ mm/d]
					//q [zt] [zn]= max_float(0,(wcs[zt][zn]-wc [zt] [zn])*dx);
				break;
				default:// �brige Kompartimente

					wcint (zn,ti);
					this->set_push_back_fluxes("q",zn,ti,							 
												Basics_math::max_(0.0,
												(this->get_states("wcs",zn,end)
												-this->get_states("wc",zn,start))
												*this->get_dx(zn)*10/ti.get_dur()));			//flux from, zn-1 to zn   [ mm/d]
					//q [zt] [zn-1]= max_float(0,(wcs [zt] [zn]-wc [zt-1] [zn])*dx);
					this->set_states("wc",zn,end,
												this->get_states("wcs",zn,end)
												-guelpprog(zn+1,ti)/this->get_dx(zn)
												*(this->mod_soiltemp->get_states("t",zn,end)>0));
					//wc [zt] [zn]= wcs [zt] [zn]-guelpprog(x,zt, zn+1)*(sote[zt][zn]>0)/dx;
					
				break;
			}//Ende switch
		} //zn		
}

double Mod_wc_fc::wcint (int layer,time_intervall &ti)
{
	//returns    the real infiltartion into the soil [cm/ time_intervall]
        timestamp start=ti.get_start();
        timestamp end=ti.get_end(); 
	double pot_infiltration=0;
	double wc;
	double dx=get_dx(layer);   //cm
	double real_infiltration;
	if (layer==0)
	{
		pot_infiltration=infprog(ti);//  [cm/Zeitintervall] da wird wc[zt][zn]abz�glich evapotransp berechnet	 für das gesamte profil	
		real_infiltration=Basics_math::min_(pot_infiltration,guelp_surf(ti));	 //[cm/intervall]
		this->set_push_back_fluxes("Real_Infiltration",ti,real_infiltration*10/dx);	 //[mm/d]
		
		wc=this->get_states("wc",layer,end);
		this->set_states("wcs",end,wc+real_infiltration/dx);
		//wcs [zt] [zn]=(wc [zt] [zn]) +min(a, guelp_surf(x,zt,zn))/dx;
	}
	else
	{
		double flux_from_above=Basics_math::max_(0.0,guelpprog(layer,ti))*(this->mod_soiltemp->get_states("t",layer,end)>0);	  //[cm/timeintervall]
		wc=this->get_states("wc",layer,end);
		this->set_states("wcs",end,wc+flux_from_above/dx); // [cm3/cm3]
		//wcs [zt] [zn]=(wc [zt] [zn]) +max_float(0, guelpprog(x,zt,zn)*(sote[zt][zn-1]>0))/dx;
	}
	return real_infiltration;
}
double Mod_wc_fc::infprog (time_intervall &ti)// Nettoinfiltration Und Abfluss bzw. Speicher (surface flow is not considered every "too much" is stored)
{		
        timestamp start=ti.get_start();
        timestamp end=ti.get_end();
		get_ETR(ti); //berechnet ETR und die Wassergehaltsverteilung nach Abzug von ETR
        double pot_infiltration=this->mod_surface_water->get_states("Surface_Storage_without_infiltration",end)*ti.get_dur()/10; //[cm/time intervall]
            
		return pot_infiltration;
	
}

double Mod_wc_fc::max_possible_evaporation_of_segment (timestamp &ts,int zn)
{
	// Returns maximum of water of one segment that can be evaporated [cm3/cm2]
	 double g=0;
         double dx;
         double wc_zt_zn=this->get_states("wc",zn,ts);
         double pwp=this->mod_pdf->get_states("FK_pf4.2",zn,ts);
         
         if (zn==0) dx=inp->horizons->at(zn);
         else dx=inp->horizons->at(zn)-inp->horizons->at(zn-1);
               
         g= Basics_math::max_( dx*(wc_zt_zn-pwp*exp(epauft*(inp->horizons->at(zn)-dx/2))) ,0.0);
	 //g= max_float( dx*(wc [zt] [zn]-1*hydro_props[x] [3][0]*exp(epauft*(zn*dx-dx/2))) ,0);
	//printf("%f  \n",g);
	return g;
}

double Mod_wc_fc::max_possible_evaporation_of_profile(timestamp &ts)
{
        // Returns maximum of water that can be evaporated [cm3/cm2]
	double g=0;
	int o=0;
	for (o=0;o<this->inp->horizons->size();o++)
	{
		g=g+max_possible_evaporation_of_segment(ts,o);	
	}
	if (g<=0)
	{
		err::err3("Mod_wc_fc::max_possible_evaporation_of_profile: Maximum_possible_evaporation_of_profile_is_negative.\n");g=10;
	}
	return g;
}
double Mod_wc_fc::get_ETR(time_intervall &ti)
{
	timestamp start=ti.get_start();
        timestamp end=ti.get_end();
        int o=0;
	double max_poss_evap=max_possible_evaporation_of_profile(start); //[cm3/cm2]
        double ETP=(this->mod_etp->get_fluxes("ETP_Evaporation",ti)+this->mod_etp->get_fluxes("ETP_Transpiration",ti))*ti.get_dur()/10;//from mm/d to cm/intervall
        double dx;
        
    double  ETR=min(max_poss_evap,ETP); //Einheit ist cm
	for (o=0;o<this->inp->horizons->size();o++)
	{
            if(o==0) dx=inp->horizons->at(o);
            else dx= inp->horizons->at(o)-inp->horizons->at(o-1);
            
            this->set_states("wc",o,end,
                    this->get_states("wc",o,start)-(ETR * max_possible_evaporation_of_segment(start,o))/max_poss_evap/dx);
		
            //wc[zt][o]= wc[zt-1][o]-(et_r[zt]*maxd (x,zt-1,o)/maxsum (x,zt-1)/dx);//Einheit cm/kompartiment
		
	}
	return ETR;
}

double Mod_wc_fc::guelp_surf (time_intervall &ti)
{
    //returns the amount of water which can be uptaken by the soil [cm]
        timestamp start=ti.get_start();
        timestamp end=ti.get_end();
	double K1=0,K2=0, dh=0,h1=0,h2=0,q=0;
        double dz=this->inp->horizons->at(0)/2;
     // Pressure computation in    
        h1=0; //atmospheric pressure is 0							   //[cm]
        h2= h_after_Williams(this->get_states("wc",0,end),0,start);   //[cm]
        dh=max(0.0,h2-h1);
	 // conductivity calculation
	
        K2=this->mod_pdf->get_states("KS",0,start)/pow((double)(h2+1),(double)1.5);	 //[cm/d]
	
	
        
	q=(K1+K2)/2*(1+dh/dz)*ti.get_dur();//[cm/per time_intervall= c,3/cm2/time_intervall]
	q=min(q,dz*2*(this->mod_pdf->get_states("Porosity",0,start)- this->get_states("wc",0,end)));
	//printf ("k, dh %i %f  %f\n",zn, (K1+K2)/2, dh);
	return q;
}

double  Mod_wc_fc::guelpprog (int layer,time_intervall &ti)// Aufnahmepotential
{
		// returns the water flow from the layer above to the recent layer [cm/ time intervall]
        timestamp start=ti.get_start();
        timestamp end=ti.get_end();    
	double K1=0,K2=0;
        double h1=0,h2=0,dh=0;
        double dz=10;
        double q=0,q1,q2;
        double dx1,dx2;
        double porosity1=this->mod_pdf->get_states("Porosity",layer-1,start);
        double porosity2;
        double fc1=this->mod_pdf->get_states("FC",layer-1,start);
        double fc2;
	if (layer==this->inp->horizons->size())
	{	   //case for calculation of fluxes from last layer to subsurface
			
                if (layer-1<=0)
                {
					err::err2("Mod_wc_fc::guelpprog: There is just one layer defined");
                    dx1=this->inp->horizons->at(layer-1);
                }
                else
                {
                    dx1=this->inp->horizons->at(layer-1)-this->inp->horizons->at(layer-2);
                }
                dx2=dx1;   
                dz=dx1;
          //pressure computation      
				h1= h_after_Williams(this->get_states("wcs",layer-1,end),layer-1,start);
                h2= h_after_Williams(fc1,layer-1,start);
                dh=max(0.0 ,h2-h1);
         //conductivity computation       
				K1=this->mod_pdf->get_states("KS",layer-1,start)/pow((float)(h1+1),(float)1.5);
				K2=K1;
                
		 //flow
				q=(K1+K2)/2*(1+dh/dz)*ti.get_dur();//[cm/per intervall]
                
				q1=dx1*(this->get_states("wcs",layer-1,end)-fc1);
                
				q2=dx2*(porosity1- fc1);
				q=Basics_math::min_(q,q1);
				q=Basics_math::min_(q,q2);
				q=Basics_math::max_(0.0,q);
		
	}
	else
	{
				fc2=this->mod_pdf->get_states("FC",layer,start);
				porosity2=this->mod_pdf->get_states("Porosity",layer,start);
				
                if (layer-1<=0)
                {
					err::err2("Mod_wc_fc::guelpprog: There are just two or less layer defined");
                    dx1=this->inp->horizons->at(layer-1);
                }
                else
                {
                    dx1=this->inp->horizons->at(layer-1)-this->inp->horizons->at(layer-2);
                }
                dx2=this->inp->horizons->at(layer)-this->inp->horizons->at(layer-1);
                dz=dx2/2+dx1/2;
             //Pressure height calculation   
				h1= h_after_Williams(this->get_states("wcs",layer-1,end),layer-1,start);
				h2= h_after_Williams(this->get_states("wc",layer,end),layer,start);
				dh=Basics_math::max_(0.0,h2-h1);
			//conductivity
				K1=this->mod_pdf->get_states("KS",layer-1,start)/pow((float)(h1+1),(float)1.5);
				K2=this->mod_pdf->get_states("KS",layer,start)/pow((float)(h2+1),(float)1.5);
		   //flow
				q=(K1+K2)/2*(1+dh/dz)*ti.get_dur();//[cm/time_intervall]
                
				q1=dx1*(this->get_states("wcs",layer-1,end)-fc1);
                
				q2=dx2*(porosity2- this->get_states("wc",layer,end));
                
				
				q=Basics_math::min_(q,q1);
				q=Basics_math::min_(q,q2);
				q=Basics_math::max_(0.0,q);
		//printf ("k, dh %i %f  %f %f\n",zn, q1,q2,q);
	}
	return q;
}



double Mod_wc_fc::h_after_Williams(double wc,int layer, timestamp &ts )
{
    double a1=this->mod_pdf->get_states("Williams_a1",layer,ts);
    double b1=this->mod_pdf->get_states("Williams_b1",layer,ts);
    
    return exp(a1+b1*log(wc));
}